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Abstract 

The aim of this paper is to discuss the appropriate modelling of in- and outflow bound¬ 
ary conditions for nonlinear drift-diffusion models for the transport of particles including 
size exclusion and their effect on the behaviour of solutions. We use a derivation from a 
microscopic asymmetric exclusion process and its extension to particles entering or leav¬ 
ing on the boundaries. This leads to specific Robin-type boundary conditions for inflow 
and outflow, respectively. For the stationary equation we prove the existence of solutions 
in a suitable setup. Moreover, we investigate the flow characteristics for small diffusion, 
which yields the occurence of a maximal current phase in addition to well-known one-sided 
boundary layer effects for linear drift-diffusion problems. In a one-dimensional setup we 
provide rigorous estimates in terms of e, which confirm three different phases. Finally, 
we derive a numerical approach to solve the problem also in multiple dimensions. This 
provides further insight and allows for the investigation of more complicated geometric 
setups. 


1 Introduction 

Transport phenomena of crowded particles and their mathematical modelling have received 
considerable attention recently, driven by a variety of important applications in biology and 
social sciences, e.g. transport of ions and macromolecules through channels and nanopores 
(cf. [4, 13, 12, 15, 21]), cargo transport by molecular motors on microtubuli (cf. [8, 27, 34]), 
collective cell migration (cf. [25, 33, 36, 14]), tumour growth (cf. [37, 24]) or dynamics of 
human pedestrians (cf. [9, 22, 35]). Such applications naturally lead to questions related to 
boundary (or interface) conditions restricting in- or outflow of particles, and the resulting 
characteristics of flow. In ion channels the characteristics are relations between bath concen¬ 
trations (boundary values) and current, in pedestrian motion one is interested in flow and 
evacuation properties depending on exit doors, and the movement of cargo between micro¬ 
tubuli respectively delivery to the desired site act as as similar boundary conditions. A variety 
of computational investigations of such issues have been performed, partly with additional 
complications such as electrostatic interactions, chemotaxis, or herding. Such simulations 
can give hints on the flow behaviour, but it becomes difficult to understand the causes and 
asymptotic regimes for certain effects. Therefore we investigate in detail a canonical simple 
model with in- and outflow boundary conditions in this paper. To this end, we assume that 
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the boundary dQ of our domain is subdivided into three parts: Inflow T, outflow E and insu¬ 
lating \ (r U E) with T n E = 0. Then for x G D C M n , t > 0 and p = p(x , f) we consider 
the equation 


d t p + V-j = 0, j = -DVp + p(l - p)u, 
with boundary conditions 


( 1 . 1 ) 


-j. n = a{l-p), on T, (1.2) 

j-n = Pp, on E, (1.3) 

j ■ n = 0, on dQ \ (T U E), (1.4) 

where u : M n —> M n is a given velocity field. The model we use is derived from the paradigmatic 
asymmetric exclusion process (cf. [7]), with appropriate modifications to include realistic br¬ 
and outflow boundaries. In a simple one-dimensional setup, this model was investigated 
recently in [38] including stochastic entrance and exit conditions, exhibiting three different 
phases of behaviour. We will take a continuum limit of that model and verify that these three 
phases are still present under the same conditions on parameters and demonstrate how the 
model generalizes to multiple dimensions and multiple species going in potentially different 
directions. The (formal) continuum limit naturally leads to the case of nonlinear convection 
dominating the diffusion, hence the limit of diffusion tending to zero is natural, and indeed 
the appearing boundary layers are separating the different phases. 

Let us mention that the study of continuum limits of microscopic particle models with 
size exclusion effects is a very timely research topic. The majority of the rigorous analysis 
is however carried out for closed systems, i.e. under no-flux boundary conditions or on the 
whole space, where such systems possess a gradient flow structure that can be well exploited 
either with transport metrics (cf. [1, 6, 28]) or with entropy dissipation techniques (cf. [5, 3]). 
The case of non-closed systems has been studied at the continuum level mainly for Dirichlet 
boundary conditions, where at least the modelling is obvious. In the case of general inflow and 
outflow conditions the modelling of boundary conditions needs to be adapted to the specific 
approach for deriving continuum equations (cf. [16]), which seems to have been ignored by 
most authors in the past. Moreover, the case of non-equilibrium boundary conditions poses 
additional challenges on the analysis, in particular existence proofs for stationary solutions 
cannot be carried out anymore by explicit computations or energy minimization arguments. 
Nonetheless, some arguments can still benefit from the underlying gradient flow structure in 
the energy, in particular a transformation to dual variables (also called entropy variables) is 
quite benefitial for existence proofs, since it yields maximum principles that do not hold for 
the original variables (cf. [3, 4, 23]). In this paper we will use similar ideas and extend them 
from Dirichlet to inflow and outflow boundary conditions. 

This paper is organized as follows: In Section 2 we present the model for several species 
and give more details about the nonlinear boundary conditions. In section 3 we present 
existence proofs for a single species. We seperately treat the cases when the velocity field is 
either a given divergence free vector field or the gradient of a potential function. In Section 4 
we investigate the behaviour for a small diffusion coefficient and compare this to the results 
presented in [38]. Finally, in Section 5, we introduce a discontinuous Galerkin scheme and 
present examples in one and two spatial dimensions. 
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2 Modelling 


Crowding models based on (totally) asymmetric exclusion processes as well as their mean-field 
continuum limits have gained strong attention recently (cf. e.g. [32, 31] and the references 
above). The main paradigm is to model jumps of particles on a discrete lattice with jump 
probabilities consisting of unoriented parts (diffusion) and oriented drifts (transport). The 
exclusion is incorporated by avoiding jumps to an occupied cell. Using standard continuum 
limits (rescaling time and space to have grid sizes and typical waiting times converge to zero) 
as well as simple mean-field closure assumptions, which can also be made rigorous (cf. [18]), 
one obtains partial differential equations of the form 


d t pi + V • (ji) = 0, ji = -Di(p 0 Vpi - piVpo) + PiPoUi, in ff x (0, T), (2.1) 


where x G 0 C R n , t > 0, pi = Pi(x,t ) is the density of the i-th species of particles with 
diffusion coefficient Di and velocity field m : R n — > R n m, i = 1 The free-space 

density po is given by 

M 

po = i _ Pi- (2-2) 

i=i 

In the previously well-investigated case of a potential field Ui = —DiS/Vi for some V) : R n R 
(cf. [3]), the system can be recast in gradient form 

dtPi = V • ( Dip 0 piV(d Pi E[pi ,... ,pm})), (2.3) 


with the entropy functional 

, • • •, Pm] = / 

Jii 

The above differential equations have been studied in detail with potential fields and no-flux 
boundary conditions, when the system is indeed a gradient flow and stationary solutions can 
be characterised as minimisers of the entropy at fixed mass (cf. [3]). In many practical 
applications different boundary conditions and non-zero flow is of fundamental importance 
however. In [4] the case of mixed no-flux and Dirichlet boundary conditions has been studied 
in a model for charged particles coupled with the Poisson equation. Here we want to focus on 
in- and outflow boundaries, as recently also used in one-dimensional stochastic models [38]. 



M 

£ 

i= 1 


(.Pi log Pi - PiVi) + Po log Po dx. 


(2.4) 


2.1 Inflow Boundary Conditions 

We assume that particles of the i-th species enter the domain fl on a subregion T,; c dVL with 
rate a* > 0. Without exclusion principle, this would simply mean in the continuum that 
the normal flux equals ttj. Modelling volume exclusion in the discrete setting means that 
the particle can only enter a grid cell adjacent to the boundary if it is empty. Hence, the 
probability of entering is modified to a t po, and we deduce the boundary condition 

—ji • n = atipo on C. (2.5) 

Note the negative sign in front of the normal flux since we use the convention of a normal 
oriented outward. The boundary condition can be rewritten as 

DiPoPiV (log — ) • n = po{ai + piUi ■ n), 

V Po J 
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which clarifies the role of the normal velocity at the inflow boundary. The inflow rate can 
balance the normal velocity only if m ■ n < 0. On the other hand we will have p\ < 1 and 
thus, balancing only for a% < 1. 


2.2 Outflow Boundary Conditions 

Outflow boundaries are more straightforward to model, we assume that particles of the i-tli 
species leave the domain 0 on a subregion T, t c <90 \ Tj with rate /%. Thus, 

ji ■ n = (5iPi on £j. (2.6) 


Again we can rewrite the boundary condition in the form 


-DipopiV (log — ) ■ n = pi(Pi — poUi • n), 

V PoJ 

hence Ui ■ n > 0 is needed for balancing. Note also that we could easily include no-flux 
boundaries by simply setting (3i = 0, but we rather shall consider them explicitly, i.e. we have 


ji ■ n = 0 



(2.7) 


3 Basic Properties 

As detailed in the introduction, we shall from now on restrict our analysis to the case of 
a single active species (M = 1). If we further assume that the diffusion coefficient D is 
normalized to 1, equation (2.1) considerably simplifies to 

dtp + V • (-V/9 + p( 1 - p)u) = 0, 

with p denoting the density of a single species and supplemented with boundary conditions 
(1.2)-(1.4). We mention that with similar arguments as in [37, 24], the case M = 1 can also 
be derived from standard continuum fluid mechanical models adding a congestion constraint 
Pq = 1 — p. Due to the boundary conditions there is obviously no mass conservation in the 
system. However, there is still a natural balance condition between in- and outflow, i.e., if p 
solves (1.1) then 

dt p dx = / (3p do — apo do, (3-1) 

Jn 7s Jr 

where a and /3 denote the in- and outflow rate for p, respectively. In the stationary case the 
two integrals need to balance, which implies an interesting coupling in the balance conditions 
(via po) if M > 1. Note also that in an evacuation case, i.e. E = 0, the mass in the system 
is monotonically decreasing as it is expected. Finally, we state the following assumptions for 
later use 

Assumption 3.1. (Al) fl C M n , n = 1,2,3 with boundary <90 of class C 2 . 

(Al’) In addition to (Al), let T and £ be such that a weak solution w of the Poisson 
equation with right-hand side in L 2 (0) and Neumann data d n w = / satisfies w £ H 2 {VL) for 
any / such that 

/Is G i7 1/2 (£), /| r £ i7 1 / 2 (T), /U\ ( ruE) = 0. 
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(A2) 0 < a < 1 and 0 < /3 < 1. 

(A3) u G [IF 1,00 (f2)] n such that V • u = 0, u ■ n = — 1 on T, u ■ n = 1 on £, and u ■ n = 0 
on dfl\ (TUE). 

(A3 1 ) u = VV, VeH 1 ^). 

3.1 Existence of Stationary Solutions 

We shall present two different proofs, one for a given velocity vector field u and a second 
one where u = VE for some potential V. in which case we can employ a transformation to 
so-called ” entropy variables”. 

Theorem 3.2 (incompressible case). Let the assumptions (Al), ( A1 ’), (A2) and (A3) hold. 
Then, the equation 


V- (-Vp + p(l-p)u) = 0, 

supplemented with the boundary conditions (1.2)-(1.4) has at least one solution u G H l (Li) n 
L°°{TL) such that 

min{a, 1 — /?} < p(x) < max{«, 1-/3} (3-2) 

Proof. The proof is based on Schauder’s fixed-point theorem. We define the set 

M = G L°°(fl) n H l (Pl) | min{a, 1 — f3} < p < max{«, 1 — P}, J \Vp\ 2 dx 

with C = IMIto + 2|<?f2|. For given p G M, we define the operator S : M —> H 2 (Q) 
p to the solution of the linearized problem 


-V • Vp + (1 

— 2p)Vp ■ u = 0, 

(3.3) 

supplemented with the boundary conditions 



Vp • n = (a — 

P)(!-Ah on T, 

(3.4) 

— Vp ■ n = ((3 — 

(1 - p))p, on £, 

(3.5) 

Vp • n = 0, 

on dn\ (ru£). 

(3.6) 


Note that we linearized the boundary conditions differently on T and E, this will be crucial 
later on. Standard theory for linear elliptic equations (and our assumption on the regularity 
of the boundary), cf. [19, 26], ensures a maximum principle and existence of a weak solution, 
subsequently the existence of a solution p G H 2 (Ll) since the prerequisites of (AT) are satisfied. 
In order to apply Schauder’s fixed-point theorem, we have to prove that S is self-mapping 
from Al to Al, continuous and compact. 

Self-mapping: Equation (3.3) satisfies a maximum principle with vanishing normal 
derivative on dLl \ (r U £), and thus (by Hopf’s maximum principle) p attains its maximum 
on T U S. We have to distinguish the following cases: 

• p attains its maximum on T and thus V p ■ n > 0. Since by assumption (1 — p) >0, this 
implies a + pu ■ n > 0. As u ■ n = — 1 on T we conclude p < a. 

• If p attains its maximum on £ we conclude, since u-n = l, p<l — f3. 


< C 
that maps 


5 



If p attains its minimum on the boundary we use the same arguments to conclude a < p 
and 1 — (3 < p. Finally, the L 2 -bound on Vp follows by using the weak formulation with test 
function p and applying the bounds 0 < p < 1 and 0 < p < 1. 

Continuity: To show continuity of S we take a sequence pf. in L°°(Q) n H l (Q) such 
that pk p- We have to show that the sequence pk = S(pk) converges to some p and that 
p = S(p). Since pk £ H 2 (Q) we know that there exists p such a subsequence that ph t —*■ p 
weakly in H 2 { fl). Thus 


Vpfc, -V(j>dx 


/ p{\ — pkj)u ■ V(j)dx —> / Vp-Vcfdx 


/ p( 1 — p)u ■ V^dx, 

Jn 


i.e. p solves (3.3). The continuity of the trace operator allows us to pass to the limit in the 
boundary conditions (3.4)-(3.6) as well. The maximum principle discussed above implies that 
this solution is unique and the uniqueness of limits therefore yields p = p. 

Compactness: The compactness of the operator S follows from the fact that the em¬ 
bedding H 2 (Q) ^ L°°(n) n H 2 (Q) is compact for n < 3. This completes the proof. □ 


Next we treat the potential case u = VV, where we obtain the following 


Theorem 3.3 (potential case). Let the assumptions (Al), (A2) and (A3’) hold. Then, the 
equation 


V • (-Vp + p(l - p)VV) = 0, xeLl, (3.7) 

supplemented with the boundary conditions (1.2)-(1.4) has at least one solution u £ H l (Ll) n 
such that 0 < p < 1. 

Proof. Our proof is based on an approximation procedure, applied to the equation in entropy 
variables which are defined as the variation of the entropy functional with respect to the 
density p. Since we are in the case M = 1, the entropy functional (2.4) reduces to 

E[p\ = [ pln(p) - pV + (1 - p) ln(l - p). 

Jn 

Then, we introduce the entropy variable as 

:= d p E[p\ = log p- logpo - V. 

Using elementary calculations, we can express the original density p as 

e i’+ v 1 

p = i + e P+v ’ and Po = i + e p+v ■ 

Applying this transformation to (3.7), (1.2)—(1.4) yields the nonlinear equation 

( e ^+v \ 

-V • -= 0, (3.8) 

\(1 + e^+vf 


6 



supplemented with the boundary conditions 


e b+v 

(1 + e ^+ v ) 2 
e*+ v ' 
(1 + e b +' / )2 

e^+ v 

(1 + e ^ +vr ) 2 


Vip.n = a i+exj)+v , 

on r, 

(3.9) 

e^ +v 



Vip ■ n = B- - r—rr, 

1 + e^+ v 

on E, 

(3.10) 

Vip ■ n = 0, on <90 \(fU E). 

(3.11) 


We will now apply an approximation procedure to this equation and proceed in several steps. 

Existence for an auxiliary problem: To simplify our notation we introduce the func¬ 
tion A(ip, V ) := -— e , 2 and for 5 > 0 we consider the problem 

(l+e ^+*') 


-V • (A(ip s , V)Vip 5 ) + Sip 6 = 0. 


(3.12) 


To prove existence of (3.12), we use a fixed-point argument. For given ip 6 L 2 (fi) we define 
Ag(x) = A(ip(x), V{x)) + 5 which yields the linear equation 

-V • (A s Vip S ) + Sip 5 = 0, (3.13) 


subject to the nonlinear boundary conditions 

l 

I ut - 

1+6 

AgVip 6 ■ n = { 


-P 


1 +e^ 6+v ’ 
J s +v 


l+e^+v ’ 


on r, 
on E, 


l 0, on dQ \ (T U E). 

The corresponding weak formulation, for <p £ H 1 (12), is given by 

r / _ \ r i r p i> s +V 

0 = / ( AgS/ip 6 ■ V<£> + Sip 6 p ) d x — a -- —(pda + (3 / -— p da. 

> Jr l + e* 5 + v ,/ s i + e i>*+v r 


(3.14) 


(3.15) 


This is the Euler-Lagrange equation to the nonlinear minimisation problem for the energy 
functional 

E(i> S ) = lJ (i5|V^| 2 + <5|^| 2 ) dx-aj F(ip 6 , V) da + f3 J G(iJ> s ,V) da, 

where F and G are chosen such that 

1 e^ +v 

d^ F (ip, V ) = l + e ^ +v , dpG{iP, V) = 1 + fiV , +y - 

Note that F, G are convex, since their second derivatives are non-negative. Furthermore, due 
to 

e i>+y i 

Ag(x) = -=—;-b S = -~-b 6, 

(1 + e^ + ^) 2 2(1 + cosh(V’ + V)) 

we have that Ag(x) £ L°°(f2), uniformly with 5 < Ag(x) < <5 + 1/4. Thus E(ip 6 ) is coercive 
with respect to the H l norm and due to its convexity we conclude (cf. [17, Section 8.2, 
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theorems 2 and 3]) the existence of a unique minimiser which is by definition a 

weak solution to (3.13). Furthermore, since , V) > 0 and F('ip s , V) < oo, we infer the L 2 
a-priori estimate 

[ |V^| 2 d x<C + a [ F(?p s ,V)da < C M - (3.16) 

Jn Jr 

Here the constant Cm depends on the geometry, a, /3, and 5. This result allows us to define 
the nonlinear operator K : L 2 (fi) —> Ff 1 (fl) mapping ^ to the solution of (3.13)-(3.14). Our 
aim is to apply Schauder’s fixed point theorem in the set 

M := {</> € L 2 (n) \ m 2 LHn) < C M }. 

To this end we denote by I H i(ci)^>-L 2 (n) the compact embedding of H 1 into L 2 and define the 
operator K : M. M by K = ° AT. Since the a-priori estimate (3.16) implies 

that K is self-mapping, it remains to show its continuity. We consider a sequence p>n that 
converges to tj) in L 2 (0). This yields a sequence '</4 £ having a weak limit. Since Ag is 

uniformly bounded in L°°(f2), an application of Lebesgues Theorem yields Ag(p n ) —> Agpijj) 
in L 2 (0) and thus we can pass to the limit in the first integral of the weak formulation (3.15), 
i.e. 


/ Ag^nWi-Vtpdx^ / AsW)Vr-V<Pdx. 
In Jn 


Uniqueness of the weak solution to (3.15) (due to the convexity of the Energy E ) thus implies 
K(i/j n ) —> K(ip) in F 2 (VL). Thus Schauder’s fixed point theorem yields the existence of a 
solution to the auxiliary problem (3.13)-(3.14). 

Limit 5 —> 0: To this end, we define 


/= e 


ip A +V 


i _|_ e r+v ■ 

Then p 5 £ H 1 (12) and satisfies the equation 

V • (-V/ + p S ( 1 - /)w) - dA'ip 6 + 8ip s = 0, 
with the boundary conditions 

f -«(1 - P S ), on T, 
(-Vp S + /(I - /)VUj n = l pp 5 , on E, 


0 , 


on an\(TUE). 


Again, we consider the weak form given by 

0 = J (\ 7 p 5 - p \ l - /)vv) • Vip + 5VJ’ 5 ■ V<p + 8%!> s tp dx 

~ a J (1 — p S )<pd<j + P J p S ipda, p£H l { 12). 


(3.17) 


(3.18) 


(3.19) 


(3.20) 



Our aim is to derive a-priori estimates on p° by choosing the test function ip = . We have 

0 = [ (Vp s ■ Vip s - p s ( 1 - p 5 )VV • Vip 5 ) d x + 5 [ IVV’' 5 ! 2 + O' 5 ) 2 dx (3.21) 

./n Jn 

— a J (1 — p S )ip S da + /3 J p 5 ip s dcr. 

We estimate each term separately, noting that V'lp 6 = ~ VV". For the first term 

we have 

f ivi 2 ^ „ r V7T r V7 „<5 j „ i [ „5/1 „5mv7t^i2j_ 


f J dx - 2 [ VV-Vp 5 dx+ [ p 5 (l -p s )\VV\ 2 dx 
n p d {l-p d ) J n J n 

// (1 V)|vv| 2 d, 


>2 [ |V//| 2 dx- i [ |VF| 2 dx, 
is) 4 Jq 

where we used Cauchy’s inequality to estimate the mixed term and the fact that //(l — //) < 

1/4. For the second term we estimate 

— a J (1 — p <5 )V ;5 da = a J ^(1 — p 5 ) log -— j— + 2 p 5 — 1^ da + a J (1 — p 5 )V + 1 — 2 p 8 ^j dcr. 

The first term in this equation is a Kullback-Leibler distance and thus non-negative. As 
V € H 1 (Q), the trace theorem yields F| 9r2 € L 2 (dQ) and since, by definition 0 < p s < 1, the 
second term is bounded. For the third term of (3.21) we write 

P j^p s ^&a = P jL 5 \ogr^+2p 8 da-0 j^(-p 5 V -l + 2p 5 ) da. 

By the same arguments as above, we conclude that the first term is non-positive while the 
second one is bounded. Summarizing, we obtain 

J |V/| 2 dx< ^ J \VV\ 2 dx + aj (-{l-p 5 )V+ l-2p 5 ^ da -/3 J [~p 5 V - 1 + 2p 5 ) dcr. 

These estimates yield a a-priori bound for p s in i7 1 (ll). Due to 0 < p s < 1, this allows us 
to pass to the limit in the weak formulation (3.19). In particular we have, by passing to 
subsequences if necessary, 


/ V/x • V<pdx —> / Vp-V<pdx 

In Jn 


since V p 5 —" Vp in L 2 (D), 


f p s (l — p S )W ■ V<pdx —> f p(l — p)W ■ V<pdx since p s -» p in L 2 (D), 

Jn Jn 


S / Vip ■ V(p + ip ip dx -» 0 

Jn 


since ip s , Vip s G L 2 (fl), 


—a / (1 — p s )ip da + j3 / p 5 ipda — a / (1 — p)<pdcr + (3 / ppda since p s —>• p in L 2 (<9D). 

Jr J s Jr Jt, 
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Note in the potential case we can only conclude that the density p takes values between 
zero and one, for a stronger result depending on the in- and outflow parameters we need to 
return to the assumptions for incompressible velocity fields: 

Corollary 3.4. Let the assumptions of theorem 3.3 and additionally AV = 0, d n V = —1 
on F, d n V = 1 on E, and d n V = 0 on OLl \ (r U S) hold. Then, the solution p to (3.7) 
supplemented with (1.2)-(1.4) satisfies the bounds 

min{o:, 1 — f3} < p(x) < max{a, 1-/3} (3.22) 

Proof. Since for AV = 0 the equation (3.7) fulfills a maximum principle, the assertion follows 
by similar arguments as in the proof of Theorem 3.2. □ 

Remark 3.5. Note that testing the weak form (3.19) with the entropy variable if s yields to 
estimates analogous to those that are obtained by the entropy dissipation method in the time 
dependent case. In fact, if equation (3.17) would feature the additional term dtp (parabolic 
case), then differentiating the entropy functional with respect to time would yield 


d t E{p 5 ) = V • ((-V/h 5 + /(l - p s ) W - 5Vip s ) + Sip 6 (In p 5 - ln(l - p s ) - V) dx. 

Jo. L 


Recalling the definition = (In p d — ln(l — p d ) — V) and after integration by parts one obtains 


d t E{p 5 ) = ~f (-V/ + P 5 { 1 - p s ) W) • - 5\\hf 6 \ 2 + 6\V 

— a J (1 — p S )fi S da + (5 J p 5 da. 


dx 


(3.23) 


This means that in the entropy dissipation, we obtain the same terms as in equation (3.21). 
While in the stationary case, their sum is zero, we can conclude boundedness in the parabolic 
case by integrating (3.23) with respect to time to conclude 

E(p s ) + J^ J (-Vp s + p s (l-p s )VV^) - 5\Vip s \ 2 + 6\fi 5 \ 2 dxdt 

I f (1 — p s )fi s da dt + (3 f f p s da dt < E(pf) < C, 

Jo Jr Jo Jt, 


— a 


/o Jr 

where po denotes the initial datum. 


Remark 3.6. We finally mention that for convenience we used a diffusion coefficient equal 
to one, but all results of this section remain true for an arbitrary value D > 0 and even for 
regidar spatially varying coefficients. This is important for the following section, where study 
the natural case of a small diffusion coefficient. 


4 Asymptotic Unidirectional Flow Characteristics 

In the following we discuss in detail the flow properties of the single species model for small 
diffusion D = e <C 1, in particular for the stationary solution p £ x L°°(0) of 

V • (— eVp + p(\ — p)u) = 0, in Q (4.1) 
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with boundary conditions (1.2)-(1.4). We are interested in the asymptotic behaviour as e f 0, 
in particular the boundary layers and the asymptotic flow patterns, which we expect to be 
characterised by three different phases as in [38]: 

• An influx-limited phase with an asymptotically low density corresponding to a density 
of outgoing particles on £. 

• An outflux-limited phase with an asymptotically high density corresponding to a density 
of incoming particles on F, with a boundary layer on £ created by lower outflux rates. 

• A maximal current phase with asymptotic density ^ and boundary layers both at £ 
and r, which occurs at high in- and outflow rates. 

First we note that a direct conclusion from the maximum principle (3.2) is the non- 
appearance of a maximal current phase if 

7 , i [min{a, (1 - /?)}, max{a, (1 - /?)}]. 

In this case the maximum principle implies that the densities are bounded away from ^ 
uniformly in e. 

4.1 Characterization of Phases in the One-Dimensional Flow 

We now turn to the one-dimensional case with constant velocity, where we can use a scaling 
of space and flow such that P = [0,1], u = 1, T = {0}, and £ = {1}. This setting corre¬ 
sponds exactly to the continuum limit of the setting in [38], and we will rigorously show that 
indeed the same behaviour as in the TASEP with stochastic entrance and exit conditions - 
with transitions between the phases at exactly the same parameter values - appears for the 
continuum limit. For convenience we restate the one-dimensional version of (4.1), (1.2)-(1.4) 
as 

~ed xx p + d x (p(l - p)) = 0 in (0,1), (4.2) 

with boundary conditions 

ed x p = (1 — p)(p — a) at x = 0, (4-3) 

ed x p = p( 1 — p — fl) at x = 1. (4-4) 

A first result particular for the one-dimensional case is the uniqueness of a solution: 

Proposition 4.1. There exists exactly one weak solution p £ H 1 (P) of (4.2)-(4.4). 

Proof. Let p\ and p 2 be two solutions, then w = p\ — p 2 satisfies 

-ed xx w + ^((1 - pi - p 2 )w) = 0 


with boundary conditions 

— ed x w + (1 — p\ — p 2 )w = —aw at x = 0, 

— ed x w + (1 — pi — p 2 )w = fiw at x = 1. 
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Now let V € H 2 {[ 0,1]) be such that —edV = (1 — p\ — p 2 ) and w = e v v. Then, v is the weak 
solution of 

d x (e v d x v ) = 0 

in (0,1) with boundary conditions 

ed x v = av at x = 0, 

ed x v = —f3v at x = 1. 

Using the weak formulation of this boundary value problem with test function v implies 

[ e v v 2 dx + ae y(0) u(0) 2 + /3e y(0) n(l) 2 = 0, 

Jo 

which yields v = 0 and thus uniqueness of the solution. 


□ 


We start our analysis of the flow properties with a simple calculation relating the difference 
of p to the constant state \ to the boundary values: 

Lemma 4.2. Let p G H 1 ([0,1]) be the unique weak solution of (4.2)-(4.4). Then the estimate 


dx + /3p( 1) — — < e| 1 — a — fd\ 


(4.5) 


holds. 


Proof. Using the test function t p{x ) = x in the weak form of (4.2) and adding and subtracting 


| we find 


(td x p + (p- ^) 2 ) dx + fJp(T) ~\ = °- 


Further integrating the first term and using the a-priori bounds from the maximum principle 
for the boundary values concludes the proof. □ 

Lemma 4.2 will yield the desired asymptotic estimate if we can guarantee that (3p( 1) > 
such that the second term on the left-hand side is nonnegative. Note also that in spatial 
dimension one the flux is constant, thus we find (3p(l) = a(l — p( 0)), i.e. the above result 
could equally be formulated in terms of a respectively the inflow boundary value. To prove 
the latter under appropriate conditions is the objective of the next result: 

Theorem 4.3 (Maximal Current Phase). Let p G if 1 ([0,11) be the unique weak solution of 
(4.2)-(4.4) and let 

min {a,P} > -. 


Then the estimate 


dx < e\ 1 — a — (3\ 


(4.6) 


holds and furthermore, we have j > 1/4. 
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Proof. Using Lemma 4.2 it suffices to show (3p( 1) >1/4, which we carry out by contradiction. 
Assume p(l) = ^ — 8 with <5 > 0. Since (3 > ^ and a > ^ we conclude in particular 




and p( 0 ) = 1 - —p( 1 ) >1 - + —8 >\ + —8. 

a 4 a a 2 a 


Let H be a smooth monotone function such that 

0) = 0, H{ 1) = 1, supp(H') C (^ — 7 , ^ + 7 ) 

with 7 < min{e>, Now we choose the test function p = H(p) in the weak form of (4.2) 
again with | added and subtracted. Then we find 

[ {tH’{p)\d x p\ 2 - H\p)p(l-p)d x p) dx +Pp(l)(H(p{l)) - H(p(0)) = 0. 

Jo 


Using the nonnegativity of the first term and rewriting the second term yields 


(ft>( 1) - j)(ff Wl)) - Wp( 0)) < - £ H\p)(p - ix = F(p( 0)) - F(p(l)), 

where F satisfies F'(p ) = H'(p)(p — ^) 2 and -F(O) = 0. With the properties of H it is 
straightforward to see that 

H{p{ 1)) - H(p{ 0)) = 1, F(p( 1)) = 0, F(p( 0)) < 3 7 2 . 


Hence, we conclude 

-5 = /3p{l) -3 7 2 , 

which is a contradiction for 7 sufficiently small. Since the flux is constant, the fact that 
j > 1/4 follows immediately from (1.3). □ 


We remark that the strategy of proof of Theorem 4.3 is reminiscent of entropy solution 
concepts for conservation laws and parabolic equations (cf. [ 20 ]), where roughly speaking the 
Heaviside function applied to p — c for arbitrary constant c multiplied with a nonnegative 
smooth function is used as a test function to define entropy inequalities. The function H in 
the above proof will indeed approximate the Heaviside function of p — ^ as 7 tends to zero. 

In the inflow- and outflow-limited case the analysis is easier, an estimate like in Lemma 
4.2 suffices: 

Theorem 4.4 (Inflow- and Outflow Limited Phases). Let p 6 1L 1 ([0,1]) be the unique weak 
solution of (4.2)-(4.4) and let 

max{a,/3} < -. 

Then for a < (3 the estimate 

[ \p-a\ dx < e 1 a —— (4.7) 

Jo P~ a 

holds, while for a > (3 

[ \p ~ 1 + P\ dx < e- — a -^-- (4.8) 

Jo a - P 
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Proof. Note that the maximum principle implies a < p < 1 — /5 in any of the two cases. 
We only detail the case ^ > (3 > a, as the other one is analogous. Using the test function 
(f(x) = 1 — x in the weak formulation and some rewriting we have 

0 = [ {ed x p + p 2 - p) dx + a{ 1 - p(0)). 

Jo 


With some rearranging and the bounds on p we have 

1 f 1 

(1 — p)(p — a ) dx < e(l — (3 — a) + a p dx — ap{ 1). 

Jo 



Since p > a we have 


(P - «) 



|p — a| dx < e(l — (5 


a). 


□ 


Summing up, we have shown exactly the same behaviour for our the continuum model as 
[38] for the discrete TASEP. 


4.2 Explicit Solutions in One Spatial Dimension 

In this section, we briefly discuss explicit solutions to the one-dimensional equations (4.2)— 
(4.4). This derivation is mostly based on basic calculus and the details are presented in the 
appendix. However, this approach allows us to clarify the role of the parameter e with respect 
to the phase diagram. In particular, we can show that for £ > 0, maximal current can occur 
for values of a, (3 that are strictly smaller than 1/2. Since in one space dimension, the flux 
j is constant, we can set j = 1/4 and integrate (4.2) and obtain the first order ordinary 
differential equation 

-d x p + p(l - p) = 


which is also known as “logistic equation with harvesting” in the context of population dy¬ 
namics, cf. [2, 10]. Solving this equation subject to the boundary conditions elementary 
calculation shows that on of the following conditions on a and (3 have to hold in order to 
obtain a continuous solution: 


1 1 + 2e 

2 4e + 1 

1 1 + 2e 

2 4£ + 1 


< a < 

<P< 


1 

2 

1 

2 


and 

and 


1 4ck£ + 2a — 1 

2 8 ae + 2a — 2e — 1 ’ 

1 4/3e + 2a - 1 

2 8{3e + 2/3 — 2£ — 1 ’ 


(4.9) 

(4.10) 


Interestingly, maximal flux is achieved for values of a, f3 < 1/2 which is in contrast to the 
discrete model, [38]. To illustrate this, we depicted the changes in the phase diagram for 
different values of e in figure 1. In section 5.2 we present numerical results based on a 
discretization of (4.2)-(4.4) that confirm this observation. 
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0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 

Figure 1: In this phase diagram, the solid lines separate the regions of j > 1/4 (maximal flow) 
and j < 1/4 for the values e = 0.01 (blue), e = 0.1 (red) and e = 0.5 (green). The dashed 
lines correspond to the discrete case of [38] where maximal flux is achieved for a, f3 > 1/2, 
only. The inlet shows a magnification of the area around the point (a,/3) = (1/2,1/2). 

5 Numerical Solution 

In this section we will describe the numerical method that we used and present some examples 
in one and two space dimensions. Our implementation is based on the discontinuous finite 
element method which is well-suited for convection dominated problems, see [ 11 ] and the 
references therein. We will not give any details regarding error estimates and the convergence 
of our algorithm which remains future work. 

5.1 Setting and Discontinuous Galerkin Scheme 

Let us recall some well-known notations and definitions, cf. [11]. We start by dividing our 
domain into elements which are triangles in two space dimensions and intervals in ID. For 
simplicity, we shall only discuss the two-dimensional case from now on. We cover the domain 
D C M 2 by a finite collection of triangles which we denote by Th, where h refers to the diameter 
of the largest triangle. Furthermore, we denote by F the mesh faces which are characterised 
by one of the following two conditions: 

1. Either, there are distinct triangles T\ and T 2 such that F = dT\ - F is a interface, 

2. or, there is T G Th such that F = dT n <912 - F is a boundary face. 

We denote by T l h the set of all interfaces, the boundary faces and by Fh the union of 
these two sets. Furthermore, hf is the normal vector of a facet, pointing outward. On Th we 
introduce the broken polynomial space 

Vh = {v € L\n) : VT G %, v\ T € V\T) }, 
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where V l (T) denotes polynomials of degree one on T. For a scalar function v, smooth enough 
for the expression v\ F for all F G T to make sense, we define interface averages and jumps 
in the following way 

{MM*) : = v \t 1 ( x ) + v \To ( x ))' f° r a - e - x *= Fi (average), 

[n]i?(.x) := v\ T (x) — v\ T2 (. x ), for a.e. x G F, (jump). 

With these definitions at hand, we can state our discontinuous Galerkin scheme. Starting 
from the weak formulation of a linearized version of (1.1) we consider 


: / Vp\7(j)dx+ / p{l — p)uV(j)dx + a / p(f)ds + (3 / pj>ds = a / J>ds, 4>h G If 1 (II), 

Jn Jn Jr Jr: Jr 


= :a(p,4>;p) 


=:a F {p,<t>) 


■ -f(o) 


(5.1) 


with p G H l {Q) n L°°(Q) given. In order to obtain a discrete solution ph G 14 we define the 
bilinear form 


a(p h , <t>h I p) = a swip (p h , <j) h ) + a upw (p h , <j> h ; p), 


with a symmetric weighted interior penalty method for the diffusion given by 

a swip (p h , 4> h ) = [ eVhPh ■ V h 4>hdx - V' e [ {{{’V h Ph}} ■ n F \4> h \ + {phj{{'Vh4>h}} ■ n F ) da 
Jn Jf 


FeF h 


+ ^2 Vj— [ \phWhjda 
^ h F Jf 


F£F h 

and a upwind scheme for the advection part 

a upw ( Ph , (f>h) = [ ~Ph(( 1 - Ph)u ■ Vhffih) dx + V] f ((1 - p h )u ■ n F {{ph}}{{4>h}} da 
J n FeF^ J F 


+ ^K 1 - Ph) u ■ n F \ f WI4Jda, 

FeFi F 


again with ph G 14 given. The local length scale h F is dehned as h F = \{h Fl +hx 2 ), where T\ 
and T 2 are the two triangles adjacent to face F. In order to obtain a solution to the original 
nonlinear problem (1.1) we employ the following semi-implicit iteration scheme: For given 
find u^ +1 G 14 s.t. 

« +1 ! </>/i) + r(a(u” +1 ,</> h ;^) + a F (n)( +1 ,( ? !>/ l )) = (ujj, <f> h ) + f((f) h ), \/J) h G V h , (5.2) 


with a relaxation parameter r > 0. Thus in each step one has to solve the following system 
of linear equations 

(M + rA)u n+1 = (. Mu n + r/), 

where u denotes the vector of coefficient of u in the linear finite element basis, A is the matrix 
corresponding to the bilinear form (a + a F ), and M denotes the mass matrix. The vector / 
stems from the term f(4>h) on the r.h.s. of (5.2) with u n being the solution of the previous 
step. In all experiments below we chose uq = 1/2 and r = 0.01. Note that this scheme can be 
interpreted as a semi-implicit time discretization of the parabolic version of (1.1) with time 
step size r. 
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Figure 2: Some results of the ID code for e = 0.1 (blue), e = 0.05 (black), and e = 0.01 
(red).Top left: a = 0.2, fl = 0.4, Top right: a — 0.4, ft = 0.2, Bottom left: a = 0.6, ft = 0.7, 
Bottom right: a = 0.5, /3 = 0.5 

5.2 Results in one spatial dimension 

In one space dimension, we used MATLAB to implement the scheme described above. We 
will present several examples in the following and consider the domain D = [0,1] discretized 
by n = 200 elements. 

Different phases 

First we present some examples to illustrate the occurrence of the three different phases 
(namely influx limited, outflux limited and maximal current) that are analysed in section 
4 (and also in [38]). We performed simulations for e = 0.1, 0.01, 0.001. For a and fl we 
chose the values 0.2, 0.4, 0.6 and 0.4, 0.2, 0.7, respectively. The numerical results confirm the 
predicted occurrence of three phases and the results are shown in figure 2. 

Maximal current for a<l/2or/3<l/2 

Here we present numerical evidence for the results of section 4.2, namely the occurrence of 
the maximal flow phase for a, fl < 1/2. We chose e = 0.01 which yields 2 ^ 2€ 1 +1 ^ = 0.4902. 
We chose a = 0.4912 and by (4.9), the corresponding ft is 0.603773585. To illustrate the case 
ft < 1/2 we simply interchange the roles of a and ft. Both results are depicted in Figure 3 
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Figure 3: For e = 0.01, the resulting density (left) and flux (right) is depicted for the values 
a = 0.4912, (3 = 0.6043 (top) and a = 0.6043, (3 = 0.4912 (bottom). The maximal flux 
j = 1/4 is observed in both cases. 


and confirm the results from 4.2. To further explore the behaviour of the flux, we used the 
discontinuous Galerkin scheme introduced above to numerically produce a phase diagram by 
sampling the values of a, f3 from 0 to 1 with a stepsize of 0.01. For e = 0.1 we compared the 
contour line j = 1/4 with (4.9) or (4.10), respectively, see figure 4. 

5.3 Results in two spatial dimensions 

In two spatial dimensions, we used the software package FeniCS, [30, 29] to implement the 
scheme described in section (5.1). We present several examples using the domain sketched in 
figure 5, i.e. a corridor of length 2 and height 1 with two entrances and two exists on each side. 
The upper entrace and exit are located at 0.65 < y < 0.85, the lower ones at 0.15 < y < 0.35. 
For each entrance, we have a different inflow rates Oj, i = 1,2, while we have /?*, * = 1,2 for 
the exits. 

Maximum principle 

In all examples in this section, we use a velocity field given as the gradient of some potential. 
From theorem 3.3 and corollary 3.4 we know that for general potentials V we only have 
0 < p < 1 while for V satisfying the assumptions of corollary 3.4 we have that min{a, 1 — (3} < 
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Figure 4: For e = 0.1, the left picture shows a phase diagram generated using the discontinuous 
Galerkin method described above for values of a, f3 = 0,0.001,0.002,..., 1. On the right side, 
contour lines for several values of the flux j are depicted. The line for j = 1/4 (blue line) is 
compared to the analytical results of section 4.2 (red circles). Compare also with figure 1. 


p(x ) < max{a, 1 — f3}. To illustrate this, let V m be given as the solution to the equation 

—AV m = 0 in 11, 
d n V m = -1 on T, 
d n V m = 1 on S, 
d n V m = 0 on <912 \ (r U £), 

with the normalization condition fg^u da(x ) = C. On the discrete level, we use a mixed 
method to discretize this equation in order to ensure the condition V • W m is fulfilled exactly 
on the discrete level. The normalization constrained is achieved by setting an arbitrary 
boundary node to zero. The resulting velocity field u m = W m is depicted in Figure 6 . 
Alternatively we chose Vfax ) = x which yields ui = VV) = (1,0)*. In our first example we 
then chose a\ = 0.2, = 0.4, fa = 0.4 and fa = 0.2 and apply our scheme with V = V) and 

V = V m , respectively. The results shown in figure 7 produce the expected behaviour, namely 
the maximum principle (3.2) only occurs for V = V m . Furthermore, the results show that 
the asymmetric in- and outflow rates indeed some of the “particles” entering at «2 move over 
two the exit at fa. This indicates that the model may be able to also predict lane formation 
in the case with more than one active species (i.e. (2.1) with M > 1), see also [35]. 

High densities and obstables 

In a second example, we explored the situation of maximal flow by using the values a\ = 0.6, 
«2 = 0.9, fa = 0.9 and fa = 0.6 and V = V/. Note that in both cases we observe on the parts 
between the in- and outflow boundaries that the maximum principle of theorem 3.2 does not 
hold since • n / 0 on the no-flux boundary. The results are shown in figure 8 . Since this 
example shows that high densities can occur between the two exits, we modify the domain 
by adding an round obstacle in front of the doors as shown in figure 9. This is motivated by 
results from models for human crowd motions where in some situations, an obstacle in front 
of the exits can improve the situation. Indeed, our results show that the densities between 
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Figure 5: Sketch of the geometry for the two-dimensional simulations. We consider a corridor 
of length 2 and height 1 and with two entrances and to exits, solated on the left and right 
boundary, respectively. 



Figure 6: The vector field 'S7V m . The non homogeneous Neumann boundary conditions yield 
a velocity field that transports density away from the entrances and towards the exits thus 
facilitating the tranport. 


the two exits decreases, however at the price of a large density in front of the obstacle itself. 
Furthermore, the transition from high to low densities observed in figure 7 is shifted towards 
the entrances. 

6 Summary &; Outlook 

In this paper we analyzed a model for crowded transport with a single active species. We 
started by giving some details about the modelling then proceeding with two existence proofs 
in the stationary case. Next we analysed the flow characteristic of our model in the case of 
small diffusion. In one space dimension, we were able to recover three different phases that 
were already observed in the stochastic model [38]. Further investigation showed however 
that the continous model can produce fluxes that exceed the value j = 1/4 which do not 
occur on the discrete level. We concluded by presenting some numerical examples in one and 
two spatial dimensions. 
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Figure 7: Simulation with rates au = 0.2, a .2 = 0.4, (3\ = 0.4 and @2 = 0.2. Above the results 
for V = Vi are shown and even though the velocity field points in x-direction only, density is 
transported towards the larger exit. For V = V m (below) one clearly sees that the maximum 
principle is satisfied while still some density is transported to the lower exit. 
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Figure 8: Simulation with V = Vi and rates a\ = 0.6, «2 = 0.9, f3\ = 0.9 and /?2 = 0.6, i.e. in 
the regime of maximal flow. 



Figure 9: Introducing an obstacle dramatically decreases the density in front of the two 
exists, however at the cost of a slightly increased density in front of the obstacle. The rates 
are a\ = 0.2, «2 = 0.4, f3\ = 0.4 and /?2 = 0.2. 
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Our analysis and especially the numerical examples in two spatial dimensions suggest that 
interesting phenomena can occur when dealing with more than one active species. Since each 
species has its own in- and outflow rate, it is not clear whether one would again observe differ¬ 
ent phases, clearly separated by certain values of these parameters. Also, to prove existence 
in the case M > 1 becomes much more involved. Regarding the numerical discretization, a 
scheme that uses the reformulated problem in entropy variables might be an alternative to 
the direct approach used here. 
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Appendix 


In this appendix we will detail the calculation of explicit solutions to the one-dimensional 
equation (4.2). Since the flux j is constant in ID, we can integrate this equation to obtain 
the ordinary differential equation 

-d x p + p{l - p) = j (6.1) 

supplemented with the boundary conditions j = a(l — p(0)) and j = (3p( 1). We shall 
separately discuss the cases of constant density, maximal flux and the general case j < 1/4: 

1. p = const: Working with (6.1), the problem reduces to an overdetermined algebraic 
system of equations, namely 


j = a ( 1 - p), 

j = Pp, 

j = p(l- p). 


This is solvable if and only if a + /3 = 1 and we obtain p = a = 1 — (3 and j = a(l — a) = 
(3(1 — (3). In particular, the maximal flux j = 1/4 is achieved for a = (3 = 1/2, only. 


2. j = 1/4: Note that (6.1) with j = 1/4 is also known as “logistic equation with harvest¬ 
ing” in the context of population dynamics, cf. [2, 10]. In this case, (6.1) becomes 

2 


which yields 


-ed x p - ( p - ^ ) =0, 


/ . 1 £ 

p(x) = X + ——, 
2 x T c 


( 6 . 2 ) 


with the constant c to be determined by the boundary conditions. We have 


i = a(l-p(0)) = aQ-^ 

7 = «1 ) = f>(\ + 


1 + c 2 


Cl = 


=> C 2 = 


4 ae 
2a - 1 
4f3e 


2/3-1 


- 1. 


In order to obtain a single, continous solution we have to ensure the two conditions 


ci = c 2 and either ci = c 2 > 0 or ci = c 2 < —1. 


Very elementary but rather tedious calculations show that this amounts to the following 
two conditions on a, (3 and e 


11 + 2 £ 
2 4e + 1 
11 + 2e 
2 4e + 1 


< a < 

< (3 < 


1 

2 

1 

2 


and 

and 


1 4ae + 2a — 1 

2 8ae + 2a — 2e — 1 

1 4(3e + 2a - 1 

2 8/3s + 2(3 - 2e - 1 


for c > 0, 
for c < — 1, 


with c := ci = c 2 . As already discussed in section 4.2, this in interesting since solution 
with maximal flux occur for values of cc, (3 < 1/2 which is in contrast to the discrete 
model, [38]. See also 5.2 for some numerical examples that confirm this observation. 


26 



3. j ^ 1/4: In this case the explicit solutions is given by 


p(x) 


1 yj 4j - 1 ( 1 y/4j -l(x + c) 

-tan- 

2 2 \2 e 


(6.3) 


The boundary conditions reduce to the following nonlinear algebraic system 


J = 2 + 


1 , y[4^\ 


tan 1 - 


1 + c) 


J=P\- 2 - 


1 


tan 1 - 


1 \J 4j - lc 


Solving these two conditions for the constant c, we finally end up with a non-linear 
equation determining j given by 


arctan 


2 j - a 




+ arctan 


( 2 j-/3 


, a V 4/ - 1/ 2e 

This equation can be solved for example by applying Newton’s method. 


Note that using these calculation, one can basically calculate the solution to (4.2) (with u = 1) 
explicitly for j = 1 /4 or p = const and for other values of j by solving a system of non-linear 
equations. 
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